Transport on coupled spatial networks 
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Transport processes on spatial networks are representative of a broad class of real world systems 
which, rather than being independent, are typically interdependent. We propose a measure of utility 
to capture key features that arise when such systems are coupled together. The coupling is defined 
in a way that is not solely topological, relying on both the distribution of sources and sinks, and 
the method of route assignment. Using a toy model, we explore relevant cases by simulation. For 
certain parameter values, a picture emerges of two regimes. The first occurs when the flows go from 
many sources to a small number of sinks. In this case, network utility is largest when the coupling 
is at its maximum and the average shortest path is minimized. The second regime arises when 
many sources correspond to many sinks. Here, the optimal coupling no longer corresponds to the 
minimum average shortest path, as the congestion of traffic must also be taken into account. More 
generally, results indicate that coupled spatial systems can give rise to behavior that relies subtly 
on the interplay between the coupling and randomness in the source-sink distribution. 

PACS numbers: 89.75.Fb, 05.40.-a, 64.60.aq 



Systems that can be represented as a group of inter- 
acting networks are found everywhere in modern life 
[3]. From so-called smart power grids — which couple 
electrical distribution networks with ICT control net- 
works [4 — to interactions between other types of critical 
infrastructure networks, such as food and water supply, 
transport, fuel, and financial transactions. Recent the- 
oretical studies on the subject have generated a great 
deal of interest by demonstrating that coupling two or 
more networks together can lead to system-wide be- 
haviour which differs fundamentally from the behaviour 
of each individual network [5HT4] . These studies describe 
essentially cascade-like processes where typically either 
inverse-percolation [6j [9j [11] or sandpile methods [10] are 
used (or variants thereof). In the former case, the robust- 
ness of the system is characterized by size of the remain- 
ing giant connected component, whilst for the latter, it 
is the size of the largest sand-cascade. In both cases, the 
quantity of interest is directly related to the topology of 
the network and does not permit any consideration of 
dynamical processes which may take place on the net- 
work. Furthermore, robustness against cascade failures 
is not the only consideration for those affected by such 
real world systems. 

One broad class of processes that occur on a network 
are general transport processes, or flows [15]. Whether 
flows of people, fluids, or electrical currents, these sys- 
tems can be characterized by specifying the topology of 
the underlying network, a source-sink distribution, and 
a dynamic (Fig. [T]). Where, to avoid confusion, we only 
imagine dynamical processes that converge to a steady 
state — resulting in a stationary distribution of flows over 
the network. Unfortunately, the methods of analysis 
mentioned above do not capture many of the typical fea- 
tures one might expect here. For example, it is easy to 
imagine a simple source-sink distribution that allows the 
network to be split into two distinct components such 
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FIG. 1: (Color online) A system made of two coupled net- 
works where the nodes of network 2 form a subset of the 
nodes of network 1. Edges of network 1 are shown in black, 
edges of network 2 are shown in red (gray offline), and nodes 
in common to both networks are considered to be coupled 
(shown by dashed lines). Highlighted in green (gray offline), 
we represent a path between two nodes, the "source" i and 
the "sink" j. 

that the flows are unaffected. In this case, the size of the 
giant component may decrease but the network is still 
operating well. With this example in mind, one question 
that arises is: how should an interacting, or coupled, 
set of flow networks be characterized, and what are the 
interesting features of such systems? Whilst any system- 
wide behaviour is intimately linked with the particular 
dynamics, some understanding can be gained by investi- 
gating the properties of simple examples that are chosen 
well enough to represent some sub-class of these systems. 
In this study, we report the results of investigating such 
a toy model, and highlight the interesting features which 
we believe might be typical of many problems in this 
class. 

Most existing studies of coupled networks focus on 
variants of the random graph [6] [TDJ [11] , primarily due 
to the simplicity with which properties can be calculated. 
However, many physical networks (i.e., electrical, trans- 
portation, ICT etc.) are spatial networks and are often 
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planar [16]. For this reason, this letter focusses on cou- 
pled planar networks. Although, for completeness, the 
example of spatially embedded Erdos-Renyi random net- 
works is discussed later alongside the results for the pla- 
nar case. In addition to this, coupled networks are usu- 
ally found to be linked by a set of nodes common to 
both networks (note that this is however not a neces- 
sary limitation of the model, but simply a more realistic 
assumption for spatial networks). For example, this is 
the case for a road network coupled to a rail or subway 
network. Here, all the nodes of the road network are 
not nodes of the rail network, but conversely, all sta- 
tions are located at points which can be considered as 
nodes in the road network. Motivated by this simple 
example, we construct a first planar network as a trian- 
gulation of points in the plane. Triangulations are often 
used as a convenient way to generate planar networks 
from a given distribution of nodes. We choose the usual 
Delaunay triangulation [17], which typically avoids slim 
triangles — not seen very often in real networks due to 
their inefficiency — and which is effectively unique for a 
given set of points. We then construct a second network 
based on a random subset of the points used to construct 
the first network. Our model thus comprises individual 
networks that are each planar Delaunay triangulations, 
forming a combined network that is not necessarily pla- 
nar (see Fig. [2]) and where the nodes of the different net- 
works with the same spatial location are linked together. 
Rather than considering a dynamical system which acts 
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FIG. 2: (Color online) Each instance of the system is gener- 
ated according to the following process: | (a) | First, Ni nodes 
(here N± =30) are positioned at random within the unit circle 
and the Delaunay triangulation is produced; (b) the second 



network is then generated by drawing N2 (here iV~2 = 10) 
nodes uniformly from the existing ones (N2 < N±) and, once 
again, computing the Delaunay triangulation; (c) the com- 
bined system is no longer planar and can be represented as a 
top-down view of Fig. [I] (where zero weights assigned to the 
dotted interconnecting lines). 

to minimize a global quantity — such as electrical net- 
works, where the dissipated power is minimized — we allo- 
cate flows on the network following a basic transportation 
analogy. Here, the source-sink distribution is replaced by 
an origin-destination (OD) matrix T^ . This has the ben- 
efit that it explicitly specifies the flow between node i and 
node j. Therefore all that remains is to decide a method 
of route assignment. The obvious choice is to use the 



weighted shortest path, where the number of such paths 
between nodes i and j is denoted by cr^-. In our model, 
the weight associated with each edge is the length of that 
edge multiplied by a factor < a n < 1 , which is common 
to all edges belonging to the same network. The subscript 
n is used to label the network: n = 1 corresponds to the 
larger network and n = 2 the smaller. The idea is that 
a = oliIol\ is a single parameter that controls the rela- 
tive weight per unit distance between the two networks. 
Indeed, in order to simplify further, we impose the arti- 
ficial constraint that a < 1. This has the effect that a 
journey on the smaller network (n = 2) is favored over 
a journey of equivalent distance for the larger network 
(n = 1). We also note that since the edge weights are 
proportional to edge length and nodes are positioned at 
random, it is very unlikely that <Jij > 1. 

Previous studies of interacting networks use the term 
coupling to describe how well two networks are linked. 
Typically, this is a purely topological definition i.e., the 
fraction of nodes from one network which link to an- 
other [5] , or the probability that a particular node has an 
edge which connects both networks [10]. For transport 
processes, a better measure of interaction must include 
details of how the flows are distributed. For the system 
outlined above we define the coupling as 
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where <r^ oupled is the number of shortest paths between 
nodes i and j, which include edges from both networks, 
and where the entries of the origin-destination matrix Tij 
are normalized i.e., . = 1. It is clear from Eq. (llj) 
that A G [0, 1] is just the fraction of travellers that use 
both networks. Such a definition is dependent on the 
method by which the flows are allocated and not just 
on the system topology. Indeed, for a given allocation 
method and network topology, there is usually a max- 
imum value of A strictly less than one. In our model, 
the coupling is controlled by choosing a. By virtue of 
changing the weights associated with each network, a 
changes the (weighted) shortest path between any two 
nodes. For example, a value of a close to one indicates 
little difference between the two networks and hence, on 
average, shortest paths do not utilize both networks. By 
contrast, a low value of a (close to zero) gives rise to 
significantly lower weights on the second network and 
therefore shortest paths typically use both networks. 

With Eq. in mind, instead of investigating the like- 
lihood of catastrophic cascade failures, we consider more 
general measures of how well the system is operating. 
For example, one such measure is the average distance 
travelled 
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where is the distance travelled between nodes i and 
j. For most practical transport processes, a well de- 
signed system reduces the average distance travelled (i.e., 
water/food supply, the Internet, transportation, etc.). 
Another important quantity, which is a simple proxy 
for traffic, is the edge betweenness centrality, defined as 
x e = ^Zi^j Tij {(Tij (e) /o~ij), where subscript e is used to 
label edges, and (e) is the number of shortest paths 
between nodes i and j, which use edge e. The between- 
ness centrality allows us to introduce a second measure 
we are concerned with, the Gini coefficient G. A number 
between zero and one, G is typically used in economics 
for the purpose of describing the distribution of wealth 
within a nation. Here it is used to characterize the dispar- 
ity in the assignment of flows to the edges of a network, 
something that has been done before for transportation 
systems such as the air traffic network [18]. For exam- 
ple, if all flows were concentrated onto one edge, G would 
be one, whilst if the flows were spread evenly across all 
edges, G would be zero. We use the definition according 
to Ref. CEI 
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where subscripts p and q label edges, E is the total num- 
ber of edg es, x p is the flow assigned to edge p as defined 
earlier, and x = x p /E is the average flow on an edge. 

Unfortunately, it is impractical to consider the inter- 
play between A, d, and G, for all possible OD matrices. 
Therefore it helps to choose a specific example. We start 
with a monocentric OD matrix — i.e., all nodes travel to 
the origin — and then add noise by rewiring in the fol- 
lowing way. For each node, with probability p, choose a 
random destination, and with probability 1 — p, choose 
the origin (see Fig. [3|. 
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FIG. 3: (Color online) Representations of OD matrices where 
each arrow corresponds to an e ntry in and which relate s to 
the set of points in Fig. [5] |(a)| A monocentric OD matrix. |(b)| 
A monocentric OD matrix randomly rewired with probability 
V 0.5. 



corresponds to an average over fifty instances of the OD 
matrix for each of fifty instances of the coupled network 
geometry. We find that, as the coupling increases, the 





0.2 0.4 0.6 0.8 1.0 

(b) 



FIG. 4: (Color online) Simulation results for the average 
shortest path and the Gini coefficient (iVi = 100, N2 = 20, 
and p values: (purple dots), 0.2 (blue squares), 0.4 (green 
diamonds), 0.6 (orange triangles), and 0.8 (red inverted tri- 
angles)). When the coupling increases, the average shortest 
path decreases and the Gini coefficient can increase for large 
enough disorder. 



length of the average shortest path decreases (Fig. 4(a) ). 
This is straightforward to understand since the increased 
coupling is simply a result of reducing a. Furthermore it 
is clear that increasing randomness in the origin destina- 
tion matrix increases the length of the average shortest 
path by an almost constant value, irrespective of the cou- 
pling. By contrast, the behaviour of the Gini coefficient 



at different couplings (Fig. 4(b) ) is less easily explained. 
Consider instead Fig. |5j Here, each colormap shows the 
distribution of flows resulting from many instances of the 
system. 




FIG. 5: (Color online) Colormaps showing normalized edge 
flows — plotted at the midpoint of each edge — over many in- 
stances of the system. Colors are assigned from lightest 
to darkest, starting with white (for zero flow) and moving 
through yellow, orange and red for higher values of flow, until 
reaching black (maximum flow). Each Subfigure corresponds 



to the following para meter values: 
p 0.2, a 0.1; 



p = 0.8, a 0.9; |(d) 



0.2, a 
• = 0.8, , 
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The set of numbers iV~i, N 2l p, and a, now define an en- 
semble of systems that are statistically equivalent (with 
respect to A, J, and G). We proceed by calculating the 
quantities (A), (d), and (G) for different values of p and a, 
where angle brackets (...) represent an ensemble average. 
The results are shown in Fig. |4j where each data point 



The first two plots, Figs. 5(a) and |5(b)| were gen- 
erated from OD matrices rewired with low probability 
(p = 0.2) i.e., almost monocentric. The ratios of edge 
weights per unit distance between the two networks are 
a = 0.9 and a = 0.1 respectively. Therefore each dia- 
gram corresponds to a point on the blue line in Fig. |4(b)| 
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For a = 0.9, there is minimal coupling between the net- 
works and a high concentration of flows are seen around 
the origin. Since the flows are disproportionately clus- 
tered, this configuration is described by a high Gini co- 
efficient. By contrast, for a = 0.1, the difference in the 
edge weights means that it can be beneficial to first move 
away from the origin in order to switch to the 'fast' (low 
a) network. We therefore see a broader distribution of 
flows with small areas of high concentration around cou- 
pled nodes. The emergence of these hotspots away from 
the center also corresponds to a high Gini coefficient — 
and therefore the blue line in Fig. |4(b)| is relatively flat. 
Figs. |5(c)] and |5(d)] correspond to the red line of Fig. |4(b)| 
generated from OD matrices rewired with high probabil- 
ity (p = 0.8). We observe that even for a close to one, 
the localization of flows is less than for p = 0.2 — resulting 
in a lower Gini coefficient. As a is decreased, the second 
network becomes more favourable and coupling hotspots 
can be seen once again — resulting in a high Gini coeffi- 



cient and a positive gradient for the red line of Fig. 4(b) 



This result points to the general idea that randomness in 
the source-sink distribution leads to local congestion and 
more generally to a higher sensitivity to coupling. 

Heuristically, one might consider this a simple model 
of a two mode transportation system in the low density 
regime — i.e., where the effects of congestion do not af- 
fect route choice. We imagine a road network coupled 
to a rail network where users select the quickest route 
to their destination. That is, the shortest path actually 
represents the quickest path. This implies that the scale 
factors a roa d and a ra n must have units of time divided by 
distance, so we assume that they represent the inverse of 
the average speed associated with each mode. The result 
is that decreasing (increasing) the ratio of these factors, 
a, has the effect of increasing (decreasing) the relative 
speed of rail above road — and hence the coupling. In 
this picture, the Gini coefficient can now be thought of 
as a measure of road use. A low value indicates that the 
system uses all roads to a similar extent, whilst a high 
value indicates that only a handful of roads carry all the 
traffic. With this analogy in mind, it is natural to com- 
bine the effects observed above into a single measure. We 
assert that it is likely that a designer or administrator of 
a real system would wish to simultaneously reduce the 
average travel time and minimize the disparity in road 
utilization. To serve this purpose, we define a 'utility' 
function F = (d) + /i(G), where it is immediately appar- 
ent from Fig. [4] that, for certain values of /i, the function 
F will have a minimum. That is, a non-trivial (i.e., non- 
maximal) optimum A will emerge. Figure 6(a)| shows 
that, whether a non-trivial optimum coupling exists de- 
pends on the origin-destination matrix. For OD matrices 
rewired with a high probability, increasing the speed of 
the rail network reduces the road utilization as flows be- 
come concentrated around nodes where it is possible to 
change modes. Dependent on the value of /i, the effect 



of reduced utilization can outweigh the increased jour- 
ney time, leading to a minimum in F. Monocentric OD 
matrices, by contrast, have inherently inefficient road uti- 
lization when applied to planar triangulations, regardless 
of the speed of the rail network. Therefore no minimum 
is observed, and hence no (non-trivial) optimum A. More 
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FIG. 6: Existence of an optimal coupling: (a) Simulation 
results for \i — 10, N\ — 100, AT 2 = 20, and p values: (purple 
dots), and 0.8 (red inverted triangles) (only two values of p 
are shown to ensure the lines of best-fit can be seen clearly) . 
[(b)] Minima of quadratic best-fit curves for different values of 
p. We obtain p* ~ 0.32 using a straight-line approximation 
between the two closest points, above and below, A* = 1. 
(The error bars shown are those of the closest data point to 
the minimum of the best-fit curve). 



systematically, we plot in Fig. 6(b) the minima A*(p) of 
quadratic best-fit curves obtained from Fig. 6(a) t each 



corresponding to a different value of p. Defining , the 
value of p for which A*(p*) = 1, it is possible to catego- 
rize the system into one of two regimes. We observe that: 
if p < p* , then the optimal coupling is trivially the maxi- 
mum; otherwise if p > p* , a non-trivial optimal coupling 
exists. 

Finally, we note that similar effects can be observed 
on other non-planar networks. As mentioned earlier, one 
example is Erdos-Renyi random networks generated for a 
set of nodes with random locations. Whilst the observed 
behaviour is qualitatively the same, the results are much 
less pronounced due to the absence of spatial structure 
and the spatial localization of centrality. 

In conclusion, the model is characterised by two com- 
peting forces — the desire to move all flows onto the most 
efficient network, whilst also ensuring that congestion 
does not arise around the nodes which connect both net- 
works. We observe that the optimisation of such a sys- 
tem can be sensitive to randomness introduced in the 
origin-destination matrix. The broader interpretation of 
our work is that, spatial, space- filling, networks such as 
transportation networks or the electricity grid, may be in- 
herently fragile to certain changes in supply and demand, 
such as the transition from centralized power generation 
to decentralized prosumers [20]. This behaviour is not 
captured by the existing literature, and demonstrates an 
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alternative view of transport processes on interacting net- 
works. Indeed, whilst the assumptions made have a con- 
venient interpretation in terms of bimodal transportation 
systems, we expect that the results hold for a broader 
class of systems and welcome work in this area. 
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